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Characterization  of  solar  collectors  is  based  on  experimental  techniques  next  to  validation  of  associ¬ 
ated  models.  Both  techniques  may  be  adopted  assuming  different  complexities.  In  this  work,  a  general 
methodology  to  validate  a  collector  model,  with  undetermined  associated  complexity,  is  presented.  It 
serves  to  characterize  the  device  by  means  of  critical  coefficients,  such  as  the  film  (convection)  transfer 
coefficient,  plate  absortance  or  emmitance.  The  first  step  consists  of  identifying  those  significant  param¬ 
eters  that  match  the  selected  model  with  the  experimental  data,  via  nonlinear  optimization  techniques, 
applied  to  steady  state  conditions.  Second,  new  correlations  must  be  adopted,  in  those  terms  where  it  is 
necessary  (i.e.  film  coefficient  equations).  Finally,  the  overall  model  must  be  checked  in  transient  regime. 
To  illustrate  the  technique,  a  tailor-made  prototype  flat  plate  solar  collector  has  been  analyzed.  An  inter¬ 
mediate  complex  collector  model  has  been  proposed  (2D  finite-difference  method).  Both  steady  and 
transient  states  were  analyzed  under  different  operating  conditions.  Parameter  identification  is  based  on 
Newton’s  method  optimization.  For  parameter  approximation,  exponential  regression  functions  through 
multivariate  analysis  of  variance  is  proposed  among  many  other  alternatives.  Results  depicted  a  robust¬ 
ness  of  the  overall  proposed  method  as  starting  point  to  optimize  models  applied  to  solar  collectors. 
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1.  Introduction  and  previous  work 

Processes  of  industrialization  and  economic  development 
require  important  energy  inputs  [1].  Fuels  are  the  world’s  main 
energy  resource  and  are  considered  the  center  of  energy  demands. 
However,  reserves  of  fossil  fuels  are  limited  and  their  large-scale 
use  is  associated  with  environmental  deterioration  [2].  These  facts 
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Nomenclature 

a/  constant  factor  of  the  regression  function  T 

bui  exponent  (constant  value)  of  the  regression  function 

T 

Bi  Biot  number 

Cf  fluid  specific  heat  (J/(kg  K)) 

cp  specific  heat  of  the  absorber  plate  (J/(kg  K)) 

Dj  inner  diameter  of  the  risers  (m) 

e  plate  thickness  (m) 

eq  2D  equivalent  thickness  of  the  pipe-weld-plate  (m) 

er  pipe  thickness  (m) 

fi  exponential  regression  function  T 

Gr  Grashof  number 

h  final  estimated  film  coefficient  of  the  inner  wall 

water  (W/(m2  K)) 

h0  initial  film  coefficient  of  the  inner  wall-water 

(W/(m2  K)) 

Ij  global  incident  radiation  over  tilted  surface  (W/m2 ) 

ke  extinction  coefficient  of  the  covering  glass  (m-1 ) 

kf  pipe  conductance  (W/(m  K)) 

kp  conductance  of  the  absorber  material  (W/(m  K)) 

kh  proportionality  factor  applied  to  h 

kq  proportionality  parameter  associated  to  the  volu¬ 

metric  equivalent  heat  transfer  coefficient  Uc 
liuc  proportionality  factor  applied  to  Uc 
kiT  conduction  heat  transfer  coefficient  of  the  insula¬ 
tion  material  (W/(mK)) 

Kc  global  heat  transmission  coefficient  plate-fluid 

(W/(m2  K)) 

rrif  fluid  flow  rate  into  the  risers  (kg/s) 

rhj  fluid  flow  rate  in  lower  and  upper  pipes  (kg/s) 

n2  refraction  index  of  the  covering  glass 

Nu  Nusselt  number 

Pr  Prandtl  number 

q'[  unitary  surface  heat  radiation  on  the  absorber 

(W/m2) 

q\ ,  5  internal  heat  generation  in  the  absorber  due  to  solar 
radiation  (W/m3) 

q  2  unitary  surface  heat  loss  in  the  absorber  plate 

(W/m2) 

q2  loss  of  the  internal  heat  generation  in  the  absorber 

(W/m3) 

q  3  unitary  surface  heat  transferred  to  the  working  fluid 

(W/m2) 

q3  loss  of  the  internal  heat  generation  in  the  absorber 

due  to  heat  transmission  to  water  (W/m3) 
qv  internal  heat  generation  in  the  absorber  (W/m3 ) 

q'u  net  heat  released  in  the  refrigerant  water  per  length 
unit  (W/m) 

Vj  inner  radius  of  the  riser  (m) 

re  outer  radius  of  the  riser  (m) 

r'i  inner  radius  of  lower  and  upper  pipe  (m) 

Re  Reynolds  number 

s  collector  tilt  (rad) 

SPn  known  conditions  vector  at  time  ‘n’ 

t  time  (s) 

T  temperature  (K) 

Ta  cuidado  si  es  T0  (Eq.  (2.b)) 

T0  room  temperature  (K) 

Tb  plate  temperature  in  the  bulb  thermometer  (K) 

Tfi  refrigerant  fluid  inlet  temperature  (K) 

Tf0  refrigerant  fluid  outlet  temperature  (K) 

Tf  fluid  mean  temperature  in  a  differential  element  (K) 


Tf*  fluid  temperature  into  the  upper  and  lower  horizon¬ 

tal  pipes(K) 

Tp  material  temperature  of  the  upper  and  lower  hori¬ 

zontal  pipes  (K) 

Tm  mean  top  surface  temperature  (K) 

Tw  mean  temperature  in  the  inner  side  of  the  pipe  in  a 
differential  element  (K) 

UPn  Unknown  parameters  vector  at  a  stationary  state  at 
time  ‘n’ 

UL  volumetric  heat  transfer  loss  coefficient  of  the  col¬ 

lector  (W/(m3  K)) 

U/  global  heat  loss  coefficient  of  the  collector 

(W/(m2  K)) 

Ub  heat  loss  behind  the  collector  (W/(m2  K)) 

Ut  top  heat  loss  (W/(m2  K)) 

Uc  volumetric  equivalent  heat  transfer  coefficient  of 
the  control  volume  pipe-weld-plate  (W/(m3  K)) 

Xui  standardized  input  parameter  ‘u’  of  a  regression 
function  for  output  parameter  T 
x  transversal  direction  distance  in  the  collector  (m) 

y  longitudinal  direction  distance  in  the  collector  (m) 

a  thermal  diffusivity  of  the  absorber  plate  material 

(copper)  (m2/s) 

olq  directional  absorbance  of  the  cooper  black  plate 

A  increment 

sg  glass  emmitance 

sp  plate  emmitance 

p  absorber  density  (kg/m3) 

Pf  refrigerant  density  (kg/m3) 

ip  latitude  (rad) 

co  timing  angle  (rad) 

8  declination  (rad) 

y  azimutal  angle  (rad) 

6  incident  angle  of  beam  irradiance  over  the  tilted  sur¬ 

face  (rad) 

&  objective  function  to  minimize  by  the  identification 

process 

r  auxiliary  expressions  for  Eqs.  (7.a),  (7.b),  (9.a)  and 

(9-b) 

Sub-index 

i  x  direction 

j  y  direction 

m  number  of  nodes  in  x  direction 

w  number  of  input  parameters  in  a  regression  function 

p  number  of  riser 

r  node  number 

Super-index 

k  time  in  transient  regime  analysis 

n  interval  of  time 


have  encouraged  growth  in  the  use  of  renewable  energy  resources 
worldwide  [3].  Solar  energy  is  considered  one  of  the  main  promis¬ 
ing  alternative  sources  of  energy  to  replace  the  dependency  on 
other  fossil  energy  resources  [4,5].  Solar  water  heating  systems 
are  very  common  systems,  extensively  used  in  many  countries 
with  high  solar  radiation  potential,  such  as  Mediterranean  coun¬ 
tries  [6].  They  are  often  viable  to  replace  fossil  fuels  used  for  many 
home  applications  [7].  Conventional  flat  plate  solar  collectors  with 
a  metal  absorber  plate  and  covers  are  used  to  transform  solar  energy 
into  heat  [8].  Optimization  of  solar  heating  systems  provides  both 
running  and  design  parameters  that  ensure  maximal  collected  heat. 


F.  Cruz-Peragon  et  al.  /  Renewable  and  Sustainable  Energy  Reviews  16(2012)  1709-1720 


1711 


Main  parameters  to  be  considered  are  related  with  both  geometry 
and  materials  properties  [9-16].  An  optimal  design  must  ensure 
minimal  working  fluid  pressure  loss,  besides  a  nearly  constant 
temperature  on  each  transversal  section  of  the  plate.  Otherwise, 
temperature  between  pipes  will  rise,  thus  increasing  radiation  heat 
loss.  To  carry  out  the  optimization  of  the  process,  a  system  sim¬ 
ulation,  including  model  validation  and  analysis  under  different 
working  conditions  are  needed.  In  this  sense,  the  first  step  consists 
of  determining  the  parameters  associated  to  heat  collection  and 
transmission.  Subsequently,  a  mathematical  model  considering  the 
following  premises  must  be  established: 

a)  Parameters  associated  to  the  physical  properties  of  the  system. 
According  to  the  literature,  main  material  properties  that  influ¬ 
ence  the  system  are  well  known  [17-19].  Conventional  absorber 
plates  usually  exhibit  either  parallel  or  serial  hydraulic  con¬ 
figurations.  These  configurations  (including  “Z”  configuration) 
besides  circulating  regimes  showing  small  Reynolds  number 
(Re)  lead  to  non-uniform  flow  distribution.  On  the  other  hand, 
heat  transfer  from  wall  to  working  fluid  leads  to  non-uniform 
transversal  temperature  distribution  on  the  plate.  Thus,  disre¬ 
gard  of  some  critical  parameters  related  to  both  fluid  flow  and 
heat  transfer  by  convection  makes  difficult  the  optimal  design  of 
the  system.  Although  primary  surface  pressure  loss  can  be  easily 
estimated  [20,21],  little  research  about  coefficients  of  pressure 
loss  in  junctions  and  bifurcations  in  collector  pipes  has  been 
done.  In  fact,  previous  works  only  consider  parameters  depen¬ 
dent  on  Re  at  different  flow  ratios  [22,23].  Flow  ratio  is  defined 
as  the  relation  between  the  discharge  in  each  riser  (being  riser 
no.  1  the  nearest  to  the  inlet  port)  and  the  total  discharge  of  the 
collector.  In  this  sense,  there  is  an  extensive  literature  devoted 
to  semi-empirical  correlations  for  film  coefficients  approxima¬ 
tion  [24-26] ;  however,  correlations  for  pressure  loss  injunctions 
and  bifurcations  are  still  needed  [27,28].  On  the  other  hand, 
extremely  small  Re  means  axial  flow  velocity  is  distorted.  This 
feature  makes  difficult  to  understand  and  approximate  the  pro¬ 
cess  [29].  Also,  flow  distribution  seems  to  be  strongly  influenced 
by  the  Grashof  number  (Gr)  and  the  pipe  leans  [30,31  ]. 

b)  Simulation  complexity.  Collector  characteristics  are  usually  well 
known  and  needed  for  model  design.  However,  system  char¬ 
acterization  may  increase  the  complexity  of  the  model.  This 
complexity  can  be  described  under  both  stationary  and  transient 
regimes  through  computational  fluid  dynamics  (CFD)  methods 
[19,32-35].  The  target  is  to  reach  a  compromise  between  the 
complexity  of  the  selected  model  and  the  accuracy  of  its  results, 
which  should  adequately  simulate  the  system  behavior. 

Single-dimensional  methods  are  able  to  simulate  with  high  accu¬ 
racy  the  majority  of  the  situations.  Nevertheless,  to  carry  out  an 
optimal  design,  two-dimensional  methods  based  on  finite  differ¬ 
ence  or  control  volume  techniques  must  be  considered,  so  a  good 
fitting  can  be  achieved  [36,37].  In  this  sense,  to  determine  the  tem¬ 
perature  distribution  over  the  plate,  an  adequate  simulation  of 
the  heating  system  entails  an  associated  mathematical  difficulty, 
considering  both  characterization  and  CFD  techniques.  In  fact,  in 
case  experimental  data  show  high  temperature  difference  between 
upper  tubes,  previous  models  are  useless  [36,38,39]. 

For  these  reasons,  validation  of  a  defined  numerical  or  math¬ 
ematical  model  involves  data  acquisition  and  identification  of 
parameters  that  show  abnormal  behavior  according  to  the  liter¬ 
ature.  Eventually,  it  can  force  to  re-define  some  of  the  equations 
associated  to  the  model. 

Thus,  the  objective  of  this  work  is  to  propose  a  complete  pro¬ 
cedure  to  optimize  the  equations  associated  to  a  collector  model. 
In  the  first  step,  the  identification  of  the  parameters  associated  to 
materials  and  heat  transfer  provides  the  system  validation  over 


the  time.  Then,  previous  values  will  help  to  provide  new  corre¬ 
lation  functions  associated  to  the  collector  behavior.  In  sum,  a  first 
approach  to  optimize  a  two-dimensional  finite  difference  model 
associated  to  a  flat  plate  solar  collector  with  low  flow  rate  is  pro¬ 
posed. 

Identification  of  a  thermo-physical  process  involves  the  deter¬ 
mination  of  its  parameters  by  solving  an  inverse  heat  transfer 
problem  (IHTP).  Inverse  problems  consist  in  the  determination 
of  the  initial  parameters  when  solutions  are  known.  This  kind  of 
analysis  has  been  proposed  by  many  researchers  during  the  past 
years  [40,41  ].  Also,  it  has  been  proved  that  solutions  to  IHTP  are 
usually  unique,  though  estimations  are  not  always  numerically 
stable  [42].  In  this  work,  an  inverse  problem  (identification  pro¬ 
cedure)  is  performed  by  means  of  iterations  over  a  direct  problem 
(system  model)  using  a  derivative  dependent  method.  The  identifi¬ 
cation  is  based  on  a  non-linear  optimization  technique  that  uses  an 
objective  function,  obtained  from  the  sum  of  squares  of  the  dif¬ 
ferences  between  the  model  and  the  experimental  values  [43,44]. 
This  quadratic  form  provides  the  best  performance  in  a  wide  range 
of  engineering  optimization  problems.  This  technique  applies  the 
second  order  Newton  minimization  method,  in  combination  with 
a  Taylor  expansion  series  of  truncated  to  the  third  derivative 
[43].  According  to  the  literature,  it  is  observed  that  the  solution  set 
(temperatures)  is  continuous  and  smooth,  without  drastic  changes 
[36,38,39].  The  same  assumptions  have  been  applied  to  heat  trans¬ 
fer  problems  [45].  For  these  reasons,  when  a  quasi-Newton  method 
truncated  to  the  second  derivative  is  selected,  the  computational 
efforts  are  reduced  while  keeping  their  accuracy  [43].  If  the  num¬ 
ber  of  measured  points  (in  terms  of  time-space)  is  not  lesser 
than  the  number  of  parameters  to  be  identified,  the  solution  of 
a  non-linear  steady  IHTP  is  unique.  The  simultaneous  identifica¬ 
tion  of  thermo-physical  characteristics,  internal  heat  source  and 
heat  flux  would  be  unique  if  either  before  or  during  the  numer¬ 
ical  experiment,  the  monotonicity  or  piecewise-monotonicity  of 
the  parameters  to  be  determined  is  found  to  be  likewise  [40].  In 
the  present  work,  none  of  these  premises  fulfil  the  problem  condi¬ 
tions  because  only  two  sensors  were  used.  For  this  reason,  results 
only  describe  a  general  approach.  Nevertheless,  specificity  could  be 
increased  by  adding  more  sensors  over  the  absorber  and  into  the 
risers. 

2.  Materials  and  methods 

2.1.  Equipment 

The  solar  collector  system,  as  shown  in  Fig.  1 ,  is  composed  of  a  set 
of  parallel  copper  pipes  under  a  ‘Z’  disposition.  They  are  attached  to 
a  copper  flat  plate  (absorber)  by  means  of  a  tin-silver  (6%)  welding. 
This  assembly  is  covered  with  a  matt  black  paint  coating,  operat¬ 
ing  as  a  selective  surface,  thus  making  possible  a  great  radiation 
collection.  The  casing  is  made  of  a  white-lacquered  aluminium 
profile.  Two  frames  of  this  type  are  superposed.  One  of  them  con¬ 
tains  white  glass,  while  the  other  one  contains  the  pipes-plate  unit. 
The  isolator  is  composed  of  two  thin  aluminium  sheets  separated 
by  an  expanded  polystyrene  layer.  It  incorporates  inlet  and  outlet 
water  thermometers  and  manometers,  next  to  a  bulb  thermome¬ 
ter  placed  on  the  absorber  plate.  Security  valves  are  attached  to  this 
equipment.  The  main  characteristics  and  material  properties  of  this 
collector  are  shown  in  Table  1 .  They  provide  an  accurate  simulation 
of  the  system  [46]. 

2.2.  Collector  model 

Analyses  must  be  carried  out  under  rigorous  model  condi¬ 
tions,  including  the  main  features  of  a  representative  dynamic 
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LEGEND 

I. -  Manometer  in  inlet  pipe 
3.-  Cooper  absorber  plate 
5.-  Isolators 

7.-  Aluminium  case 

9.-  Manometer  in  oulet  pipe 

II. -  Thermometer  of  glass  coveri  ng 


2.-  Thermometer  in  inlet  pipe 
4.-  Pipes  and  risers 
6.-  White  glass  covering 
8.-  Purging  valve 
10.-  Thermometer  in  outlet  pipe 
12.-  Thermometer  of  flat  plate 


Fig.  1.  Collector  under  study. 


Table  1 

Collector  properties. 


Property 

Value 

Absorber  material 

Copper 

Absorber  thickness  (mm) 

2 

Absorber  useful  surface  (m2) 

0.21 

Inner  diameter  of  the  horizontal  tubes,  R'T  (mm) 

20 

Inner  diameter  of  the  riser,  R,  (mm) 

10 

Number  of  risers 

15 

Riser  length  (mm) 

450 

Distance  between  risers  (mm) 

30 

Glazing  thickness  (mm) 

4 

Complete  useful  surface  (m2) 

0.37 

Posterior  insulation  (mm) 

10  (aluminium 
sandwich)  +  7 
(polyurethane) 

Collector  plate  efficiency  factor,  FR 

0.7 

simulation  model  and  applied  only  to  the  plate-risers  unit  [47].  In 
this  particular  case,  the  following  conditions  were  adopted: 


y  a 


Nodes  i  =  1  to  7 


Fig.  2.  Sample  network. 


depending  on  the  fluid  temperature  [49].  The  necessary  data 
are  collected  from  the  literature  [50]. 

(iii)  Heat  loss  coefficients  of  the  plate  are  considered  to  be  depen¬ 
dent  on  both  average  and  room  temperatures,  at  each  interval 
of  time. 

(iv)  Temporary  variations  of  irradiation,  room  and  water  inlet  tem¬ 
peratures  are  considered  model  inputs.  In  this  sense,  both  the 
unitary  heat  release  of  the  tilted  surface  (adjusted  with  the 
help  of  an  isotropic  model,  with  reflectance  equal  to  0.2)  and 
the  absorptance-transmittance  product  depend  on  time  [17]. 
Additionally,  time  variations  of  water  outlet  temperature  and 
bulb  temperature  (over  the  flat  plate)  are  known  but  only  used 
during  the  identification  process. 

(v)  According  to  Streeter  et  al.,  flow  distribution  into  each  pipe 
may  be  calculated  by  means  of  the  friction  loss  and  conven¬ 
tional  techniques  [20].  However,  secondary  loss  coefficients 
applied  to  collector  junctions  and  bifurcations  are  both  Re  (of 
the  main  pipe  flow)  and  flow  rate  dependent  [22].  According 
to  Wietbrecht  et  al.,  they  may  be  approximated  by  exponential 
functions  [22]. 

(vi)  Film  heat  transfer  coefficients  are  initially  approached  by 
means  of  semi-empirical  equations  [24,25,28,50,51  ].  Although 
there  are  many  particular  correlations,  those  listed  in  Incropera 
et  al.  have  been  selected  [50].  Once  the  optimization  process  is 
finished,  the  coefficients  are  empirically  evaluated  to  validate 
the  model. 


(i)  A  heat  transfer  analysis  by  finite  difference  approach  of  the 
absorber  plate  and  the  fluid  that  flows  through  the  pipes  must 
be  carried  out.  Some  simplifications  consider  the  plate  as  a  two- 
dimension  system  and  each  pipe  as  a  one-dimension  system. 
The  model  consists  of  a  set  of  non  linear  equations  for  both 
steady  and  transient  states,  which  must  be  solved  by  iterations 
[48]. 

(ii)  Conduction  coefficients  and  diffusivities  of  the  plate,  pipes 
and  cover,  as  well  as  insulating  materials  must  be  taken  into 
consideration.  On  the  other  hand,  densities,  specific  heat,  vis¬ 
cosity,  conductivity  and  Prandtl  number  of  the  circulating  fluid 
(water)  are  approximated  by  4th-order  polynomial  functions, 


To  carry  out  the  analysis,  the  absorber  surface  was  divided  into 
multiple  similar  units,  using  a  mesh,  so  it  could  be  considered  a  dis¬ 
crete  problem.  A  grid  sample  is  shown  in  Fig.  2.  According  to  this, 
the  complete  study  could  be  approximated  by  means  of  finite  differ¬ 
ence  technique  [52].  Because  pipes  distribution  along  the  absorber 
plate  is  considered  to  be  constant,  only  a  portion  of  the  plate  was 
examined.  The  plate  is  attached  to  15  identical  tubes,  symmetri¬ 
cally  distributed  regarding  they-axis,  thus  dividing  the  plate  in  30 
similar  zones.  However,  only  half  of  the  area  between  tubes  is  con¬ 
sidered,  provided  the  symmetry  in  the  geometric  design.  For  this 
reason,  only  14  similar  zones  were  considered,  as  shown  in  Fig.  2,  in 
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Fig.  3.  Flow  distribution  ratio  on  each  riser. 


addition  to  another  two  half  zones  in  both  plate  extremes,  regard¬ 
ing  x-axis.  The  boundary  line  is  considered  adiabatic,  due  to  the 
presence  of  the  isolator.  Boundary  loss  is  around  3%  of  the  net  heat 
released  [49],  so  it  was  neglected. 

Flow  distribution  along  the  risers  can  be  considered  constant 
unless  the  regime  is  laminar  [22],  which  happens  when  the  flow 
rate  is  very  low  [46,49].  In  this  work,  Re  is  smaller  than  100  along 
the  circuit  of  water,  due  to  low  flow  rate  and  non-developed  flow, 
making  viscous  forces  more  important  than  inertia  ones.  Primary 
hydraulic  surface  loss  in  laminar  flow  can  be  calculated  by  the 
Poiseuille  equation,  but  form  loss  (secondary  loss)  analysis  needs 
the  application  of  some  experimental  coefficients  [20].  With  the 
help  of  these  coefficients,  pressure  loss  and  flow  distribution  along 
the  pipes  may  be  analyzed.  Another  secondary  kind  of  coefficient 
takes  extremely  high  values,  for  Re  below  100.  So,  it  makes  pos¬ 
sible  to  calculate  form  loss  by  a  general  exponential  expression 
depending  on  Re  and  the  flow  ratio  in  junctions.  It  derives  from 
the  empirical  and  validated  functions  observed  by  Weitbrecht  et  al. 
[22]. 

The  Fourier  law,  related  to  heat  conduction  through  a  material, 
given  by  Eq.  (1),  describes  the  system  dynamics  as  follows: 

qv  dT  d2T  d2T  d2T  qv  1  dT 

Oi  AT  -f - —  —  =>•  +  Tj  y  +  Tj  T  +  7 —  —  TT- 

pCp  dr  dx 2  dy2  dz2  kp  a  dr 


with 


a 

kp 


1 

pep 


where  T  is  the  temperature  (K),  t  is  the  instantaneous  time  (s), 
kp  is  the  conduction  coefficient  of  the  absorber  (W/(mI<)),  a  is  its 
thermal  diffusivity  (m2/s),  p  is  its  density  (kg/m3),  cp  is  its  specific 
heat  value  (J/(kgI<)),  and  qv  is  the  internal  heat  generation  in  the 
absorber  (W/m3). 

In  fact,  the  dependent  variables  presented  in  Eq.  (1)  must  be 
optimized  by  a  method  that  will  be  explained  in  detail  later.  The 
flow  rate  differs  along  each  riser,  as  shown  in  Fig.  3.  The  main 
consequence  is  the  different  temperature  distribution  between  the 
sections  of  the  plate.  So,  the  developed  model  of  heat  transfer  must 
be  adapted  to  each  riser  attached  to  the  flat  plate.  For  this  reason, 
the  study  is  performed  in  14  similar  portions  of  pipes  of  the  same 
geometry  (including  two  half  zones  in  both  plate  extremes),  but 
with  different  water  flow  rate  due  to  hydraulic  unbalance  caused 
by  the  high  form  loss  coefficients  associated  to  the  laminar  regime. 

The  complete  heat  balance  considers  boundary  conditions,  that 
is,  the  input  and  output  of  energy  to/from  the  plate  at  each  time 
interval  and  per  surface  unit  (W/m2).  Firstly,  the  unitary  surface 


heat  radiation,  q"  (W/m2),  is  determined  from  the  global  incident 
radiation  over  the  tilted  absorber,  IT  (W/m2),  reduced  by  the  prod¬ 
uct  transmittance-absorptance  of  the  lra'  glass-absorber  unit,  as 
shown  in  Eq.  (2.a). 

q'i(t)  =  lTra(t)  (2.a) 

Secondly,  the  unitary  heat  loss  of  the  plate,  q"  (W/m2),  is  calcu¬ 
lated  from  the  temperature  difference  between  the  selected  point 
and  the  environment,  as  follows  (Eq.  (2.b)): 


q''(t)  =  Ul(t)(T(t)-Ta(t))  with  Ui  =  Ui{T(t),Ta(t)},  (2.b) 

The  global  heat  loss  coefficient,  Lf/  (W/(m2  K)),  is  composed  of  two 
factors:  (a)  the  loss  behind  the  collector,  Ub  (W/(m2  I<)),  due  to 
the  heat  conduction  through  the  isolator  and  the  casing,  and  (b) 
the  convection  loss  to  the  atmosphere.  Its  values  depend  on  the 
material  properties  and  wind  speed.  The  top  loss,  Ut  (W/(m2  K)), 
originates  from  solar  radiation  and  re-radiation,  system  tilt,  orien¬ 
tation  and  plate  temperature. 

Finally,  q3"  (W/m2)  (associated  to  pipes  geometry)  provides  the 
available  heat  that  passes  through  the  working  fluid,  named  use¬ 
ful  heat.  It  may  be  calculated  by  means  of  the  film  coefficient,  h 
[W/(m2  K)],  related  to  the  convection  heat  transfer  from  the  inner 
wall  of  the  pipe,  with  mean  temperature  Tw  (I<)  and  water  mean 
temperature  Tf  (K),  along  a  differential  length  Ay  (m)  in  direction 
y,  as  shown  in  Eq.  (2.c). 


( ihfCf  A  Tf  +  Aynri2  pfCf(dTf  /dt)) 


(27rn  Ay) 


(2.c) 


where  riif  is  the  fluid  flow  rate  (kg/s),  T0  is  the  room  temperature  (K), 
Cf  is  the  water  specific  heat  value  (J/(kg  K)),  pj  is  the  density  (kg/m3), 
ATf  is  the  increment  of  temperature  (K)  and  rz  is  the  internal  radius 
of  each  riser  (m).  Finally,  this  procedure  comprises  some  heat  input 
to  the  water  through  the  horizontal  manifolds  before  and  after  the 
risers.  For  their  determination,  Eq.  (3)  is  included.  This  equation 
provides  the  inlet  temperature  value  into  each  riser. 

The  very  low  water  flow  rate  and  the  small  scale  provide  a  non- 
developed  flow  into  the  pipes.  Re  is  small  enough  to  distort  the  axial 
speed  of  the  fluid  along  the  riser  by  buoyancy  forces  [29].  The  filed 
flow  is  strongly  influenced  by  Gr  and  the  elevation  angle  of  the  riser 
[30].  It  makes  impossible  to  assure  the  capability  of  the  equations 
provided  by  the  literature  for  the  determination  of  the  film  coeffi¬ 
cient,  h ,  only  by  CFD  techniques.  Provided  the  difficulty  to  carry  out 
an  initial  approach  to  achieve  an  accurate  value  of  h,  it  is  calculated 
by  means  of  an  initial  value,  h0,  obtained  from  semi-empirical  rela¬ 
tions  of  heat  transfer  theory  and  multiplied  by  a  proportionality 
term  associated  to  the  coefficient,  kht  as  shown  in  Eq.  (3).  To  assure 
coherent  results,  the  optimal  value  for  the  coefficient  is  needed. 

h{t)  =  1 0/<h /i0(t);  h0  =  h0{kf ,  q,  Nu}  (3) 

On  the  other  hand,  when  Biot  number,  Bi,  is  evaluated  ( Bi  =  he\kp , 
where  e  is  the  thickness  of  the  plate),  it  gives  a  value  much  smaller 
than  one.  This  means  the  temperature  across  the  thickness  of  the 
plate  can  be  considered  constant,  giving  the  same  results  in  both 
sides  of  the  plate,  as  found  by  other  researchers  [52].  So,  the  unit 
plate-pipes  can  be  analyzed  as  a  two-dimension  system.  Accord¬ 
ing  to  this,  the  internal  heat  generation  of  Eq.  (1)  results  from  the 
boundary  conditions  of  Eq.  (2),  adapted  to  Eq.  (4),  as  follows: 

qv(t)  =  qi(t)- q2(t)~  q3(t)  in  W/m3;  q1(f)  =  S= 

q2(t)  =  and  UL  =  ^  (4) 


As  can  be  seen,  this  simplification  leads  to  the  introduction  of 
a  volumetric  heat  transfer  loss  coefficient,  UL  (W/(m3I<)),  being 
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necessary  to  analyze  in  detail  the  useful  heat.  As  a  result  of  the 
proposed  two-dimension  problem,  the  temperature  of  the  unit 
plate-pipes  across  its  section  is  assumed  to  be  constant.  Concern¬ 
ing  the  risers,  this  assumption  is  not  accurate  because  an  inner  wall 
temperature  for  each  riser,  Tw,  is  established.  This  temperature  is 
different  from  the  mean  top  surface  temperature,  Tm ;  temperature 
set  varies  from  the  top  to  the  inner  wall.  The  weld  that  connects 
the  pipes  to  the  plate  is  also  included  in  the  study.  These  consider¬ 
ations,  besides  the  geometry  and  materials  of  each  control  volume, 
difficult  the  calculation  of  the  global  heat  transfer  coefficient  of 
the  plate-fluid,  I<c  [W/(m2  K)].  So,  the  equation  of  useful  heat  col¬ 
lected  by  the  working  fluid  must  be  recalculated.  The  simplification 
from  three  to  two  dimensions  obliges  to  accommodate  the  prob¬ 
lem.  For  this  purpose,  an  equivalent  thickness,  eq  (m),  in  addition 
to  a  volumetric  equivalent  heat  transfer  coefficient,  Uc  [W/(m3I<)], 
that  considers  the  heat  from  the  top  to  the  inner  wall,  are  evaluated. 
These  terms  must  include  a  proportional  coefficient  depending  on 
the  pipe  geometry,  plate  weld  and  position,  called  kuc.  As  happened 
to  kh  term,  the  optimal  value  that  provides  an  adequate  behavior 
of  the  collector  depends  on  time.  So  the  heat  released  to  the  fluid 
matches  Eq.  (5): 


q'3(t)  =  I<c(t)(Tm(t)-Tf(t))^q3(t)  = 


<?3  (0. 


e9  = 


(e  Ax  +  7r(r2  -  r2)) 


Ax 


;  e  =  re  -  r,;  Kc(t)  =  10f 


UC 


1 


-1 


Uc  = 


=  10' 


UC 


Ax 


Q 


1 


(e  Ax  +  7i(r2  -  r2)) 
-l 


re  re 

+  -r^ln 


rz  h(t)  kp  \rt 


(5) 


Eqs.  (4)  and  (5)  help  to  determine  the  temperature  in  particular 
points  of  the  plate.  The  one  dimensional  useful  heat  that  passes 
to  the  water,  qu'(t)  (W/m),  considering  steady  state  conditions  and 
Ay,  may  be  calculated  according  to  Eq.  (6)  as  follows: 


q'uM 


Ay 


™/cf  A7> 

Ay 


27rrz/i( t)(Tw  -  Tf ) 


2.2. 2 .  Steady  state  equations 

For  each  (i  j)  point  of  the  mesh,  and  depending  on  the  increments 
Ax  and  Ay  on  the  absorber  fin,  the  numeric  adaptation  of  the  set 
of  Eqs.  (l)-(6)  for  steady  state  analysis  follows  the  expressions: 


considering  both  limits  (r0j  and  Tm+i  correspond  to  fluid  tempera¬ 
tures  associated  to  nodes  r  =  1  and  r  =  m,  respectively). 


r,Jr  +  T0j_1-TOJ(i  +  n  =  o; 


TmJr  +  T, 


-7m-,j(l+r')  =  0; 


2jrr,  Ayh 
2ril/,rC/  ' 


(7-b) 


The  set  of  (m  +  2 )n  Eqs.  (7.a)  and  (7.b)  is  used  to  establish  the  basic 
model  in  the  plate-pipes  unit,  considering  the  unknown  tempera¬ 
tures  on  the  plate  (two-dimension)  and  the  working  fluid  passing 
through  the  tube  (one-dimension).  This  is  referred  to  the  system 
portion  described  above  (Fig.  2). 

A  non-linear  set  of  equations  that  must  be  solved  by  iteration 
is  presented.  For  this  purpose,  plate  temperature  is  considered  to 
be  room  temperature.  Then,  the  rest  of  the  parameters  such  as  UL, 
Uc,  h  and  S  will  be  calculated.  A  matrix  calculus  is  built-up,  result¬ 
ing  in  a  new  set  of  plate  temperatures.  This  process  continues  until 
temperatures  converge  to  constant  values.  It  is  very  important  to 
advise  that  this  iterative  process  belongs  to  the  optimization  pro¬ 
cess  explained  below,  considering  known  all  the  parameters  that 
define  UPn. 

To  evaluate  the  set  of  Eq.  (7)  for  each  riser,  the  inlet  water  tem¬ 
perature  must  be  estimated.  For  this  purpose,  a  light  modification 
is  adopted  and  applied  to  the  input  horizontal  manifold.  In  this 
pipe,  the  temperature  is  evaluated  for  only  one  dimension  (x-axis), 
between  the  junctions  that  connect  two  consecutive  risers.  Pro¬ 
vided  that  the  global  input  flow  temperature  and  flow  distribution 
throughout  the  unit  are  known,  the  equations  set  can  be  success¬ 
fully  solved.  It  gives  information  about  the  heat  delivered  to  the 
water  in  the  lower  pipe  (with  inner  radius  r\),  from  the  riser  p  to 
the  riser  p  + 1,  separated  Al.  It  results  in  the  following  expression: 


Tp~  i  +  Tp+ -[  —  Tp 

Ax'  Ay 


2+^  Ax' A/ 
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+  Tf,P^Ax'Ay' 
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(s  +  ult0) 


v 


2nr'  Alh 

Tf,p  r  +  T/>p_ ,  -  r/iP+1  (1  +  n  =  0;  r  =  ; 


(8.a) 

(8.b) 


where  Ax'  =  Al,  Ay/=2r/z.  The  fluid  that  flows  through  the  lower 
pipe  is  non-uniform  and  continuously  reduced,  as  shown  in  Fig. 
3.  Once  the  evaluation  of  Eq.  (7)  is  finished,  the  upper  pipe  is 
analyzed.  This  will  provide  the  outlet  global  temperature  of  the 
fluid. 

The  proposed  numerical  analysis  has  been  compared  to  other 
analytical  methods  from  the  literature,  providing  similar  results  for 
the  film  coefficient  h  and  same  qualitative  tendencies  for  pressure 
loss  [53]. 


general  form  :  TMJ  +  Ti+hj  +  TUj_x  +  TZJ+1  -  TUj 
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adiabatic  fin  boundary  form  :  2TZ+1  J  +  Ti,j- 1  +  Ti,j+ 1  -  Ti,j 
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For  a  finite  different  approach,  m  corresponds  to  the  number  of 
nodes  in  ‘x’  direction  (i  varies  from  1  to  m).  Moreover,  the  work¬ 
ing  fluid  receives  a  certain  quantity  of  heat  when  passing  from  y=j 
to  j+  2.  It  flows  into  the  riser  placed  at  each  extreme  of  the  mesh, 
as  shown  in  Fig.  2.  Thus,  an  additional  equation  for  each  one,  i.e. 
Eq.  (6),  is  needed.  A  numerical  approach  is  presented  in  Eq.  (7.b), 


2.2.2.  Transient  state  equations 

When  Eqs.  (7.a),  (7.b),  (8.a)  and  (8.b)  are  used  under  transient 
state  condition,  the  inclusion  of  settings  that  consider  the  deriva¬ 
tive  of  temperature  vs.  time  is  required.  To  evaluate  it,  an  implicit 
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method  has  been  proposed.  The  method  is  very  stable  and  does 
not  require  extremely  narrow  intervals  of  time.  In  this  case,  a  simi¬ 
lar  equations  set  depending  on  known  values  of  characteristics  and 
temperatures  of  a  previous  state,  k ,  is  introduced  [52  ].  Once  temper¬ 
atures  are  known  for  t  =  k,  they  will  be  calculated  for  t=k+ 1  using 
an  additional  coefficient.  The  general  coefficients  from  Eqs.  (7.a) 
and  (7.b)  are  changed  according  to  Eqs.  (9.a)  and  (9.b),  respectively 
(other  terms  and  equations  are  adapted  accordingly). 


Ti-h  +  Cj  +  +  C i  -  Tu'A  =  s  +  ™  -  <7 


A  = 


4  + 


UL  A  a 

—  Ax  Ay  +  r 
kp 


r  = 


Ax  Ay 
a  At 


(9. a) 


fk+l  r»  ,  rrk+\ 
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rpk+ 1 

i±l  ,j 


(i  +  r  +  rp) 


—  _  r  Tk 


Ttrj2  Ay  pf  m 
4irif  At 


for  i  =  1  and  i  =  m 


(9.b) 


2.3.  Solution  to  model  optimization  procedure 

The  proposed  procedure  consists  in  the  following  steps: 

a)  The  system  and  its  model  are  defined  analytically  and  numer¬ 
ically  for  one  and  two  dimensions  [46,54].  Initially,  the  values 
of  certain  physical  properties  of  the  collector  materials  (input 
parameters  of  the  model)  are  provided  according  to  the  lit¬ 
erature.  These  parameters  are  integrated  into  the  vector  of 
characteristics  UP  that  may  vary  for  each  instant  of  time  n 
(named  UP  ).  Initially,  the  parameters  included  in  this  vector 
are  the  following: 

-  ke :  extinction  coefficient  of  the  covering  glass  (nrr1 ) 

-  n2:  refraction  index  of  the  glass 

-  olq\  directional  absorbance  of  the  copper  black  plate,  consid¬ 
ering  an  incident  angle  of  irradiation  over  the  tilted  surface,  6 , 
depending  on  tilt,  latitude,  declination,  solar  time,  etc. 

-  sp\  plate  emissivity 
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-  £g\  glass  emissivity 

-  1<iT :  conduction  heat  transfer  coefficient  of  insulation  [W/(mI<)] 

-  kp\  conduction  heat  transfer  coefficient  of  copper  [W/(mI<)] 

-  kh :  proportionality  parameter  associated  to  the  tube  inner  wall 
film  coefficient  (convection  heat  transfer),  h0  [W/(m2  I<)]  [46]. 

-  kq :  proportionality  parameter  associated  to  the  volumetric 
equivalent  heat  transfer  coefficient,  Uc  [W/(m3  K)],  consider¬ 
ing  two-dimensional  systems  [46]. 

b)  Once  the  system  is  defined,  operating  conditions  are  tested  two 
hours  a  day.  These  conditions  include  different  days  and  collec¬ 
tor  tilts.  For  each  instant  of  time  ‘n\  the  vector  SPn  contains  the 
following  input  variables  associated  to  both  the  irradiation  and 
the  system  behavior  (i.e.  flow  rate  and  temperatures): 

-  irif :  water  flow  rate  (kg/s) 

-  s:  collector  tilt  (rad) 

-  (p:  latitude  (rad) 

-  oo:  timing  angle  (rad) 

-  <5:  declination  (rad) 

-  y:  azimutal  angle  (rad) 

-  0:  incident  angle  of  beam  irradiance  over  the  tilted  surface 
(rad) 

-  It',  global  irradiance  over  the  tilted  collector  (W/m2). 

-  Ta:  room  temperature  (K) 

-  Tfii  fluid  inlet  temperature  (K) 
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Fig.  4.  Flow  chart  of  the  identification  process  (steps  a-c). 


Fig.  5.  Temperature  distribution  on  the  absorber:  (a)  isometric  view;  (b)  plant  view. 
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Table  2 

Characteristic  angles  of  the  collector  for  some  random  measures. 


Stationary  test  no.  Collector  tilt  s  (°) 

Latitude  (p  (°) 

Solar  hour  angle  co 

(°) 

Declination  8  (°) 

Azimutal  angle  y  (°) 

Incident  angle  of  direct 
irradiation  6  (°) 

1 

38 

38 

15 

20.73 

0 

25.56 

2 

28 

38 

22.5 

21.27 

0 

24.47 

3 

48 

38 

15 

21.43 

0 

34.9 

Table  3 

Input-output  data  values  of  random  measures  from  Table  2. 

Stationary  test  Solar  incident 

Room  temperature 

Inlet  water 

Outlet  water 

Bulb  plate 

Fluid  flow  rate 

no. 

radiation  IT  (W/m2) 

To  (°C) 

temperature  7,  (°C) 

temperature  T0  (°C) 

temperature  Tp  (°C) 

mf  (kg/h) 

1 

936.8 

23.2 

31 

50.5 

57 

6.42 

2 

900.8 

19.4 

26 

39 

49.5 

7.35 

3 

862.1 

25.3 

33 

46 

50.5 

7.35 

Table  4 

Optimized  values  for  unknown  parameters  after  the  identification  process  shown  in  Tables  2  and  3. 

Test  no. 

M  mm-1)  n2 

ae 

sp 

kn  (W/(mK)) 

kp  (W/(mK)) 

Uc  (W/(m3 1<)) 

h  (W/(m2  K)) 

Initial 

0.4  1.526 

0.801 

0.039 

0.88 

0.029 

402.4 

11.8  x  103 

11 

Initial  value  for 

Initial  value  for 

test  No.  1 

test  No.  1 

1 

0.4022  1.5649 

0.7543 

0.0391 

0.8804 

0.0293 

457.7 

1 1.9  x  104 

58.84 

2 

0.4016  1.5667 

0.7405 

0.0391 

0.8802 

0.0292 

423.71 

9.1  x  104 

35.7 

3 

0.4019  1.5569 

0.7668 

0.0391 

0.8803 

0.0293 

474.38 

18.4  x  104 

69 

c)  To  undertake  the  simulation,  the  variation  of  the  input  data  of 
the  vector  SPn  vs.  time  is  required.  On  the  other  hand,  vector 
UPn,  built-up  with  the  previously  mentioned  unknown  terms 
(provided  by  the  literature),  could  lead  to  a  weak  approxima¬ 
tion  to  reality.  Thus,  the  target  of  the  application  at  this  step  is 
to  identify  the  real  values  of  those  unknown  terms,  simulating 
the  real  behavior  of  the  collector  over  the  time.  This  method 
is  used  for  steady  state  analysis.  The  objective  function  &  is 
calculated  from  the  sum  of  squares  of  differences  between  the 
model  response  and  the  empirical  system  response  values,  asso¬ 
ciated  to  both  the  outlet  temperature  of  the  working  fluid  from 
the  collector  Tf0  (K)  and  the  plate  temperature  of  a  represen¬ 
tative  point  of  the  plate,  Tb  (K),  placed  near  to  the  top  of  the 
absorber  (the  real  value  is  measured  by  a  bulb  thermometer),  as 


shown  in  Figures  1  and  2.  So,  the  objective  function  remains  as 
follows: 

l //  =  ( T^al  -  T™odel  j2  +  ( j-real  _  j model )2  (10) 

If  hydraulic  friction  loss  and  fluid  flow  rate  distribution  are  con¬ 
sidered,  then  the  objective  function  must  incorporate  quadratic 
terms  related  to  pressure  in  some  points  of  the  circuit.  This 
additional  approach  has  been  previously  evaluated  in  the  same 
system  but  working  as  a  thermosyphon,  providing  accurate 
results  [53].  The  identification  procedure  is  applied  n  times  to 
determine  the  parameters  values  of  UPn,  thus  minimizing  the 
described  function  The  identification  process  described  in  this 
step  has  been  summarized  in  Fig.  4. 
d)  Next  step  comprises  the  determination  of  SPn  components  (i.e. 
heat  release  by  incident  radiation  and  angles  related  to  solar 


Fig.  6.  Temperature  variation  (°C)  in  three  cross  sections  of  the  absorber,  transver- 
sally  to  fluid  movement  at  different  locations  of  the  flat  plate. 


FLOW  TEMPERATURE  INTO  RISERS 


Fig.  7.  Water  temperature  evolution  through  pipes. 
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True  Skew 


Fig.  8.  Comparison  between  real  and  predicted  values  for  parameters  kuc  (a)  and  Nu 
(b). 

time)  and  the  system  response  data  (i.e.  plate  temperatures) 
that  fit  the  real  values  of  the  parameters  belonging  to  UPn .  If  the 
previous  expressions  do  not  suffice  to  approximate  the  desired 
functions,  other  kind  of  equations  must  be  taken  into  consid¬ 
eration,  i.e.  a  multivariate  analysis  depending  on  the  operating 
variables  included  in  SP  [55].  Statistical  fitting  functions  for 
these  parameters  must  be  analyzed.  These  functions  depend 
on  the  operating  parameters  included  in  the  SPn  vector,  mak¬ 
ing  them  time  independent.  Finally,  the  initial  model  must  be 
redefined. 

e)  Once  the  model  is  redefined,  the  transient  state  model  is  checked 
for  each  complete  daily  test. 

3.  Results 

This  section  includes  validation  of  the  numerical  model,  param¬ 
eters  identification  procedure  and  transient  regime  evaluations. 
Separation  between  nodes  in  each  of  the  14  mentioned  regions  is 
calculated  considering  Ax  =  Ay  =  0.005  m,  resulting  in  m  =  7.  This 
means  there  are  637  nodes  (7x91)  over  each  absorber  fin  (Fig. 
2).  In  the  plate  extremes,  only  one  riser  and  half  absorber  fin  are 


Fig.  9.  Experimental  and  modelled  results  at  28°  tilt:  (a)  water  outlet  temperature 
evolution:  (b)  absorber  plate  temperature  evolution. 

considered.  So,  the  number  of  nodes  is  364  (4  x  91).  Also,  the  fluid 
temperatures  in  the  boundary  limits,  considering  the  x-axis,  have 
been  considered.  Thus,  the  set  of  equations  is  composed  of  810 
equations,  with  810  unknown  values  of  temperatures  linked  to  each 
analyzed  absorber  portion.  The  intermediate  space  between  nodes 
in  the  upper  and  lower  tubes  is  0.03  m. 

3.2.  Steady  state  parameter  identification 

Initially,  the  unknown  vector  is  defined  with  data  from  the  lit¬ 
erature.  Parameters  identification  process  of  UPn  is  carried  out 
randomly  over  the  time.  Data  belonging  to  vectors  SPn  are  depicted 
in  Tables  2  and  3,  while  Table  4  shows  the  optimization  results  of 
UPn  under  the  same  conditions. 

Fig.  5  shows  the  temperatures  distribution  over  the  plate  in  a 
particular  case  study,  while  Fig.  6  exhibits  three  different  transver¬ 
sal  sections.  In  agreement  with  the  literature,  the  variation  is  small, 
but  increases  with  the  distance  from  the  riser  [17].  The  bulb  ther¬ 
mometer  shows  same  values  for  modelled  and  real  temperatures. 
Also,  extreme  values  of  input  and  output  water  temperatures  match 
real  ones.  Fig.  7  depicts  one-dimensional  temperature  variation  of 
the  fluid  flow  into  each  riser. 

As  a  result  of  the  first  evaluation,  it  can  be  seen  that  material 
properties  show  values  close  to  those  initially  estimated.  Varia¬ 
tions  in  the  glass  covering,  isolation  and  copper  plate  properties 
demonstrate  that  values  from  the  literature  are  appropriate  to  the 
simulation  process  [54].  In  fact,  they  almost  satisfy  the  condition  of 
monotonicity  needed  to  assure  an  unique  solution  of  the  IHTP  [40]. 
For  this  reason,  these  values  are  considered  known  and  removed 
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Fig.  10.  Experimental  and  modelled  results  at  38°  tilt:  (a)  water  outlet  temperature 
evolution:  (b)  absorber  plate  temperature  evolution. 


from  the  vector  of  unknown  parameters  UPn.  Thus,  computational 
charge  can  be  reduced  significantly  in  the  next  step.  On  the  other 
hand,  parameters  associated  to  heat  transfer  experiment  important 
variations  compared  to  the  initial  estimation  (from  5  to  1 0  times  the 
original  value)  provided  by  Eqs.  (3)  and  (5)  due  to  model  simplifica¬ 
tions.  The  same  conclusion  can  be  drawn  when  analytical  functions 
are  applied  [49].  Only  two  terms  belonging  to  UPn,  namely  kuc  and 
kh,  vary  Uc  and  h  values,  respectively,  over  the  time.  Flow  regime 
with  a  very  low  Re  makes  difficult  to  calculate  the  flow  rate  into 
each  riser,  the  heat  transfer  mechanism  and  the  real  convection 
coefficient.  This  model  assumes  a  two-dimension  problem,  consid¬ 
ering  the  pipes  as  flat  surfaces.  An  accurate  study  of  the  associated 
terms  Uc  and  h,  considering  physical  and  geometric  properties,  is 
required. 

Next  step  consists  in  the  identification  of  the  data  collected  two 
hours  per  day,  during  several  days.  As  mentioned  above,  only  the 
parameters  kuc  and  kh  belong  to  the  vector  UPn  and  will  be  opti¬ 
mized.  Also,  there  is  a  relation  between  Uc  and  h,  as  shown  in  Eq. 
(5). 

3.2.  New  functions  and  model  review 

Once  the  selected  parameters  are  identified,  a  relation  between 
Uc  and  temperature  is  observed.  So,  next  step  requires  to  find  out 
a  regression  function  for  parameter  /<uc,  providing  the  equivalent 
heat  transfer  coefficient,  Uc.  With  the  help  of  a  multivariate  analysis, 
the  parameter  is  adjusted  by  an  exponential/potential  function,  f\ 
[55], 


Fig.  11.  Experimental  and  modelled  results  at  48°  tilt:  (a)  water  outlet  temperature 
evolution:  (b)  absorber  plate  temperature  evolution. 

To  approximate  each  parameter  (output  variable),  regression 
functions  with  exponential/potential  shape,  //,  are  required.  For 
each  output  function,//,  it  can  be  found  wl  different  input  variables, 
Xuh  conforming  the  vector  SPn.  The  function  is  as  follows: 

fl  =  a,Xfyxfy...Xfy;  Vue(l  ,w)^Xules  (11) 

where  aj  and  bul  are  constants  of  the  function  when  it  is  successfully 
developed. 

Initially,  it  is  assumed  that f\  depends  on  different  temperatures, 
namely  inlet  and  outlet  fluid,  room  and  plate  temperatures.  Nev¬ 
ertheless,  this  first  approach  is  not  enough  to  ensure  a  good  fitting. 
So,  other  parameter  like  solar  irradiation  over  the  tilted  surface  is 
needed.  According  to  this,  the  best  result  for/i  [showing  a  standard 
deviation  of  0.033  of  the  root  mean  square  (RMS)  error  applied  to 
the  function]  is  found  when  the  released  solar  heat  under  differ¬ 
ent  temperatures  is  taken  into  consideration,  as  shown  in  Eq.  (12) 
(where  a\  is  constant). 

kuc  =  /i  ( 0\  •>  Tp,  Tjf,  Tf0 ,  T0 ,  f )  (12) 

In  relation  with  the  film  coefficient,  h,  it  can  be  seen  that  final  data 
present  a  great  dispersion.  This  confirms  the  necessity  of  establish¬ 
ing  a  new  function  to  approximate  Nu  and  identify  the  real  film 
coefficient  h  over  the  time.  It  shows  there  is  a  non-developed  fluid 
flowing  along  the  risers,  even  presenting  internal  recirculation,  as 
indicated  by  Zghal  [29].  So,  a  deep  analysis  of  the  involved  pro¬ 
cesses,  aided  by  CFD  techniques,  is  required.  However,  this  situation 
goes  beyond  the  objectives  of  this  article.  So,  a  regression  function 
to  estimate  h  and  kuc  is  enough  for  the  proposed  study. 

In  the  present  work,  kh  helps  to  determine  the  real  value  of 
h  and  subsequently  Nu.  Then,  the  procedure  for  searching  an 
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approximation  function  f2  by  stochastic  methods  is  repeated  and 
Nu  is  calculated  by  the  input  variables  Xu2 .  Between  these  variables, 
the  non-dimensional  parameters  Re,  Pr  and  Gr  provide  important 
information  about  water  flowing  into  the  risers  and  inner  wall  tem¬ 
peratures.  Nevertheless,  the  simplification  and  the  non-developed 
fluid  flow  could  lead  to  misstatements  during  heat  transmission 
modelling.  So,  the  inclusion  of  other  parameters  is  required.  As  a 
result,  the  adequate  function  includes  the  input  parameters  shown 
in  Eq.  (13)  (where  a2  is  constant),  besides  irradiation  and  tilt.  This 
last  term  affects  the  heat  distribution  and  the  transversal  velocity 
of  the  fluid  flow.  So,  determination  of  h  into  the  risers  can  be  easily 
undertaken  using  expression  (13),  as  follows: 

h=^2RL'  Nu=f2{a2,Re,Pr,Cr,Tp,Tfl,Tfo,To,It,0y,  (13) 

Results  show  a  standard  deviation  of  the  RMS  error  of  0.11.  Fig.  8 
depicts  the  correlation  between  the  functions  fa  and  f2.  Standard 
deviations  are  very  low,  thus  indicating  results  are  in  agreement 
with  parameters  values  from  Eqs.  (12)  and  (13). 

Approaches  provided  by  Eqs.  (12)  and  (13)  allow  the  validation 
of  the  proposed  numerical  model  with  available  data.  Nevertheless, 
more  general  and  useful  correlations  should  be  considered  in  future 
works,  thus  extending  the  experiments  to  fluid  flows,  temperatures 
and  friction  losses  data  [24,25,28,50,51].  In  this  sense,  increasing 
the  number  of  sensors  is  crucial.  Also,  to  improve  the  model,  Eq. 
(3)  must  be  substituted  by  Eq.  (13),  while  kuc  must  be  determined 
using  Eq.  (12)  and  applied  to  Eq.  (5). 

3.3.  Transient  state  analyses  further  to  model  review.  Final 
validations 

Once  the  model  has  been  revised,  the  transient  state  is  ana¬ 
lyzed.  Considering  the  input  data,  distributions  are  in  agreement 
with  real  data,  as  shown  in  Figs.  5-7,  but  depending  on  both  time 
and  previous  state.  From  these  distributions,  two  data  are  selected 
and  compared  with  real  values:  outlet  water  temperature  and  bulb 
temperature  (located  on  the  plate).  Some  results  are  shown  in  Figs. 
9-11. 

Approaches  present  good  accuracy.  Although  certain  delays  are 
observed,  they  are  caused  by  the  time  interval.  This  problem  can  be 
fixed  by  increasing  the  computing  time,  but  it  does  not  guarantee 
better  approximations.  In  any  case,  the  previous  model,  in  addi¬ 
tion  to  new  approximation  functions  ensures  enough  accuracy  to 
simulate  the  real  behavior  of  this  particular  system. 

Similar  results  are  found  by  other  researchers,  considering  the 
transient  behavior  of  the  collector  [35],  as  shown  in  Figs.  9-11. 
Although  results  validate  the  proposed  model,  a  deeper  study  of  the 
thermo-hydraulic  mechanisms  may  be  considered  in  future  works. 
Nevertheless,  it  is  important  to  notice  that  the  target  of  this  study 
was  to  estimate  the  temperature  distribution  over  the  plate  and 
into  de  pipes,  and  this  objective  has  been  successfully  reached. 

4.  Conclusions 

A  general  strategy  that  helps  to  validate  a  solar  collector  model 
where  some  parameters  are  difficult  to  determine,  under  both  non¬ 
steady  and  non-uniform  states,  has  been  proposed.  In  this  case,  the 
specific  conditions  of  the  fluid  regime  and  the  geometry  of  the  col¬ 
lector  make  difficult  to  understand  the  dynamic  behavior  of  the 
system.  To  avoid  the  use  of  complex  numerical  techniques,  such 
as  CFD  or  finite  element  method  (FEM),  a  highly  accurate  model 
including  simplifications  of  the  physical  unit  (plate-pipes)  has  been 
developed.  As  an  advantage,  this  process  only  requires  a  few  sen¬ 
sors.  For  this  reason,  it  makes  possible  the  analyses  in  equipments 
without  many  sophisticated  sensors  (as  used  in  heat-pipes  stud¬ 
ies).  Temperatures  distribution  of  both  absorber  and  water  gives 


qualitatively  similar  results  to  previous  works  found  in  the  litera¬ 
ture.  Moreover,  both  steady  and  transient  regimes  can  be  validated 
with  experimental  results,  thus  demonstrating  the  accuracy  of  the 
model.  Results  also  show  that  materials  properties  are  in  close 
agreement  with  the  values  given  by  the  literature.  This  means  low 
uncertainties  if  similar  values  are  used,  thus  avoiding  the  need  of 
highly  accurate  initial  data. 

Apart  from  the  type  of  model,  both  the  identification  process  and 
the  approximation  functions  constitute  two  fundamental  steps. 
Nevertheless,  the  processes  depend  on  the  following  conditions: 

1.  The  identification  process  is  based  on  a  non-linear  optimization 
procedure  over  a  quadratic  continuous  function  that  is  a  deriva¬ 
tive  dependent  function.  In  case  the  objective  function  exhibits 
a  high  discontinuity,  the  proposed  methodology  could  be  inef¬ 
fective.  If  so,  to  optimize  the  function,  other  techniques  such  as 
Genetic  Algorithms  must  be  evaluated. 

2.  Considering  parameters  approximation  functions,  exponential 
regression  functions  have  been  developed  by  multivariate  anal¬ 
ysis  of  variance.  However,  other  functions  may  also  be  evaluated, 
i.e.  linear  functions  and  potential  functions.  Including  some 
modifications  to  the  exponential  regression  functions  provides 
another  possibility  to  take  into  consideration.  In  this  case, 
stochastic  models  cannot  be  used.  So,  the  aid  of  CFD  techniques 
is  recommended. 

The  proposed  methodology  allows  the  analysis  of  the  system  under 
a  huge  range  of  different  conditions,  such  as  materials  and  geome¬ 
try  changes.  Also,  it  makes  possible  to  find  out  the  optimum  value  of 
the  pipes  diameter  and  the  distance  between  them  to  maximize  the 
captured  energy,  among  many  other  possibilities.  A  further  analy¬ 
sis  to  optimize  the  geometry  and  dimensions  may  be  developed.  In 
this  sense,  a  future  target  could  be  to  modify  them  to  maximize  the 
collection  of  energy,  as  mentioned  in  the  introduction  of  this  paper. 
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